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The Design of Approximate Hilbert Transform 

Pairs of Wavelet Bases 

Ivan W. Selesnick, Member, IEEE 



Abstract — Several authors have demonstrated that significant 
improvements can be obtained in wavelet-based signal processing 
by utilizing a pair of wavelet transforms where the wavelets form 
a Hilbert transform pair. This paper describes design procedures, 
based on spectral factorization, for the design of pairs of dyadic 
wavelet bases where the two wavelets form an approximate Hilbert 
transform pair. Both orthogonal and biorthogonal FIR solutions 
are presented, as well as IIR solutions. In each case, the solution de- 
pends on an allpass filter having a flat delay response. The design 
procedure allows for an arbitrary number of vanishing wavelet mo- 
ments to be specified. A Matlab program for the procedure is given, 
and examples are also given to illustrate the results. 

Index Terms — Dual-tree complex wavelet transform, Hilbert 
transform, wavelet transforms. 



I. Introduction 

THIS paper describes design procedures, based on spectral 
factorization, for the design of pairs of dyadic wavelet 
bases where the two wavelets form an approximate Hilbert 
transform pair. Several authors have advocated the simul- 
taneous use of two wavelet transforms where the wavelets 
are so related. For example, Abry and Flandrin suggested 
using a Hilbert pair of wavelets for transient detection [2] and 
turbulence analysis [1]. Ozturk et al. suggested it for waveform 
encoding [14]. They are also useful for implementing complex 
and directional wavelet transforms. Freeman and Adelson 
employ the Hilbert transform in the development of steerable 
filters [5], [20]. Kingsbury's complex dual-tree DWT [8], [9] is 
based on (approximate) Hilbert pairs of wavelets. The steerable 
pyramid and the dual-tree DWT have numerous benefits, 
including improved denoising capability and the fact that they 
are both directional and nearly shift-invariant. The paper by 
Beylkin and Torresani [3] is also of related interest. 

One could start with a known wavelet and then take its Hilbert 
transform to obtain the second wavelet; however, in that case, 
the second wavelet would not be of finite support. One could 
design a finitely supported wavelet to approximate the infinitely 
supported Hilbert transform, but in this paper, we design both 
wavelets together to better utilize the degrees of freedom. 

Using the infinite product formula, it was shown in [18] that 
for two orthogonal wavelets to form a Hilbert transform pair, 
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the scaling filters should be offset by a half sample. In [18], a 
design problem was formulated for the minimal length scaling 
filters such that 1) the wavelets each have a specified number of 
vanishing moments (K), and 2) the half-sample delay approx- 
imation is flat at u) — 0 with specified degree (L). However, 
this formulation leads to nonlinear design equations, and the ex- 
amples in [18] had to be obtained using Grobner bases. In this 
paper, we describe a design procedure based on spectral factor- 
ization. It results in filters similar to those of [18]; however, the 
design algorithm is much simpler and more flexible. 

A. Preliminaries 

Let the filters /io(n), hi(n) represent a CQF pair [22]. That is 

5>(n)fc 0 (n + 2*) = «(*) = 

n 

and hi(n) = (— l) n /io(M - n), where M is an odd integer. 
Equivalently, in terms of the Z-transform, we have 

H 0 {z)H 0 (l/z) -h H 0 (-z)Ho(-l/z) = 2 

and 

H x (z) = tf 0 (-l/*). 

Let the filters go(n),gi(n) represent a second CQF pair. In this 
paper, we assume hi(n) 7 gi(n) are real-valued filters. It is conve- 
nient to write the CQF condition in terms of the autocorrelation 
functions defined as 

Ph(n) := ]P h 0 (k)h 0 (k - n) = h 0 (n) * h 0 (-n) 

k 

Pg( n ) : = ^9o(k)go{k - n) = g 0 {n) * go(-n) 

k 

or equivalently as 

P h (z) := H 0 (z)Ho(l/z) 
P g (z) := G 0 (z)Go(l/z). 

Then, ho(n) and go(n) satisfy the CQF conditions if and only 
if ph(n) and p g (n) are halfband filters: 

, . fl, 71 = 0 

and similarly for p g {n). This can be written more compactly as 
p h {2n) = S{n) and p 9 {2n) = 6(n). (1) 
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The dilation and wavelet equations give the scaling and wavelet 
functions 

4>h(*) = V5 J] Ao(n)0h(2t - n) (2) 

n 

M*) = hi(n)<t> h (2t - «)• (3) 

n 

The scaling function <j> g {t) and wavelet ip g (t) are defined simi- 
larly but with filters go(n) and 

Notation: The Z-transform of h(n) is denoted by H(z). 
The discrete-time Fourier transform of h(n) will be denoted 
by H(u), although it is an abuse of notation. The Fourier 
transform of ^(t) is denoted by V(u) = T{tp(t)}. 

B. Hilbert Transform Pairs 

In [18], it was shown that if H 0 {u) and G 0 (to) are lowpass 
CQF filters with . • ■ / 

G 0 (u;) = ^ 0 (^)e"*" j( ^ ) for |u| < tt (4) 

then the corresponding wavelets form a Hilbert transform pair 

iMO = «{iMt)}- ■ 

That is 

C. Flat-Delay Allpass Filter •* . 

The design procedure presented in this paper depends on the 
design of an allpass filter with approximately constant fractional 
delay. Several authors have addressed the design of allpass sys- 
tems that approximate a fractional delay [11], [12], [16]: The 
following- formula for the maximally flat delay allpass filter 
is adapted from Thiran's formula for the maximally flat delay - 
allpole filter [23]. The maximally flat approximation to a delay 
of r samples is given by 

A( , z" L D(l/z) 

^ = —d$T ... 



or equivalently 

A(u) « e~ jujT around u = 0. 

The coefficients d(n) in (5) can be computed very efficiently 
using the following ratio. 

rf(n + l) = (n + l) (r-L) n+1 (r + l) n 
d M ^ (r-L) n . > + l) n+1 

- • _ (L - n)(L — n - r) 

" (n + l)(n>l+T) "'- 

From this ratio, it follows that the filter d{n) can be generated 
as follows: ' ■- " . " 



where 



with 



D(z) = 1 + ^2 d ^ z ' 



n=l 



J W (r + l)» 



(5) 



where (x) n represents the rising factorial 

(x) n := (x)(x -hi).. . {x + n - 1) . 
^ — v ' 

n terms 

With this D{z), we have the approximation 
A(z) « 2~ T around z — 1 



d(o) = i 

This can be implemented in Matlab with the commands 
n = 0:L-l; 

r = [1,(1, -n). * {L - n - tau),/(n + l)/(n + 1 + tau)]; 
d = aimprod(r) . 

In our problem, we will use d(n) in (5) with r =' 1/2. For 
example, with X = 2, we have 

. d(n) = {1, 2, 1/5}, for n = 0, 1,2. . - 

With L = 3, we have ' - ' 

d{n) = {1,5,3,1/7}, "fern = 0,1,2,3. 

II. Orthogonal Solutions 

In this section, we look for pairs of orthonormal wavelets* 
where me lowpass scaling filters have me form 

ho(n) = /(n) * d[n) ■ ' . . . 

9o{n) = /(n) * d(L - n) 

where the filter d(n) will be chosen to achieve the (approximate) 
half-sample delay. The first step of the design procedure will be 
to determine the appropriate filter d(n) to achieve the desired 
relationship between h 0 (n) and g 0 (n). In terms of the transfer 
functions, we have 

. H 0 (z) = F(z)D(z) 

G 0 (z) = F(z)z~ L D(l/z). 

H 0 (z) and G 0 (z) have the common divisor F(z), which will be 
determined later. We can write 



Go(z)=H 0 {z)^ 



D{z) 



where we can recognize that the transfer function 
' , z~ L D(l/z) 
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is an allpass system |^4(u;)| = 1. Therefore 

|Gb(«)| = \H 0 {w)l \Gi{u>)\ = IffiMI 

and 

|*»| = |* fc («)|. 

If the allpass system A(z) is an approximate half-sample delay 

A(iv) w c~' JUJ / 2 around a; = 0 

or equivalently, A(z) « 2 -1 / 2 around z = 1, then the sought- 
after approximation (4) is achieved 

Go(u>) « H 0 (ij)e~ j ^ around cj = 0. 

This is achieved by taking the allpass filter d{n) in (5) with r = 
1/2. (The point z = 1 is chosen for the point of approximation 
for the lowpass filter because that is the center of the passband.) 
To obtain wavelet bases with K vanishing moments, we let 

F{z) = Q(z)(l + z- 1 ) K . 

Then 

H 0 (z)=Q(z)(l + z- l ) K D(z) (6) 
Goto = Qto(l + z- l ) K z- L D(l/z). (7) 

We now have the following design problem. Given D(z) and K, 
find Q{z) of minimal degree such that ho{n) and go(n) satisfy 
the CQF conditions (1). Note that with (6) and (7), h 0 (n) and 
go(n) have the same autocorrelation function 

P(z) :=P h (z) = P 9 {z) 

= Q(z)Q(l/z)(z + 2 + z- l ) K D(z)D{l/z). 

Similar to the way Daubechies wavelet filters are obtained, we 
can obtain Q(z) using a spectral factorization approach as in 
[22]. The procedure consists of two steps. 

1) Find r(n) of minimal length such that 

a) r(n) = r(-n), and 

b) R{z)(z + 2 + z~ l ) K D(z)D(l/z) is halfband. 
Note that r(n) of minimal length will be supported on the 
range (1 - K - L) < n < (K + L - 1). 

2) Set Q(z) to be a spectral factor of R(z) 

R(z) = Q(z)Q(l/z). (8) 

To carry out the first step, we need only solve a system of linear 
equations. Defining 

S(z) := (z + 2 + z- l ) K D{z)D{\/z) 

we can write the halfband condition as 

6(n) = [l 2](a*r)(n) 

= £s(2n-A;)r(A;). 

A: 

When written in matrix form, this calls for a square matrix of 
dimension 2(K + L) - 1, which has the form of a convolution 
(Toeplitz) matrix with every second row deleted. 



TABLE I 
Matlab Program 



function [h.g] ■ hwlet(K.L) 

7. Approximate Hilbert pair of orthogonal wavelets 
X h, g - scaling filters of length 2*(K+L) 
1 K - number of zeros at z=-l 
% L - degree of fractional delay 

n ° 0:L-1; 
tau - 1/2; 

d- cumprodCCl, (L-n) .*(L-n-tau) ./(n+1) ./(n+l+tau)3) ; 
si = binom(2*K f 0:2*K); 
s2 ■ conv(d,d(end:-l : 1)) ; 
s = conv(sl,s2) ; 
M - K+L; 

C » convmtx(s' f 2*M-l) ; 

C « C(2:2:end. :); 

b - zeros(2*M-i,i); 

b(M) - 1; 

r - (C\b)'; 

q = sfact(r) ; 

f = conv(q,binom(K,0:K)) ; 

h » conv(f ,d) ; 

g « conv(f ,d(end:-l:D) ; 



The second step assumes R(z) permits spectral factorization, 
which we have found to be true in all our examples. With Q(z) 
obtained in this way, the filters Hq{z) and Goto defined in (6) 
and (7) satisfy the CQF conditions and have the desired half- 
sample delay. Note that Q(z) is not unique. 

This design procedure yields filters ho(n) and #oto of 
(minimal) length 2(L + K). A Matlab program to imple- 
ment this design procedure is given Table I. The commands 
binom and sf act for computing binomial coefficients and 
performing spectral factorization are not currently standard 
Matlab commands. They are available from the author. 

Example 1A: With K = 4 and L = 2, the filters h 0 (n) and 
go{n) are of length 12. Fig. 1 illustrates the solution obtained 
from a mid-phase spectral factorization. The plot of the function 
l^hM+j^f? tol shows that it approximates zero for u < 0, as 
expected if tph, and tp g form a Hilbert transform pair. The coef- 
ficients are given in Table II. The Sobolev exponent (reflecting 
the smoothness of the wavelets) was found 1 to be 1.983. Note 
that the Sobolev exponent is the same for tph(t) and ip g (t). 

Example IB: We set K = 4 and L = 2 again, but this 
time, we take a minimum-phase spectral factor rather than a 
mid-phase one. The wavelets obtained in this case are illustrated 
in Fig. 2. The function \^h(^) + j^g(u)\ is exactly the same 
as in Example 1 A; using a different spectral factor in (8) does 
not change it. We will see below that the mid-phase and min- 
imum-phase solutions can lead to different results when they are 
used to implement two-dimensional (2-D) directional wavelet 
transforms. (The Sobolev exponent is the same as for Example 
1A.) 

Example 2: With K = 3 and L = 3, the filters h 0 (n) and 
go(n) are again of length 12. Fig. 3 illustrates a solution using a 
mid-phase spectral factor. It can be seen that l^to +jtyg(<*>)\ 

! The Sobolev exponents were computed using the Matlab program sobexp 
by Ojanen (http://www.math.rutgers.edu/~ojanen/). 
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Fig. 2. Example IB. Approximate Hilbert transform pair of orthonormal 
wavelet bases with N = 12, = 4,L = 2, and minimum-phase spectral 
factorization. ■ ... ■ 
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Fig. 1. Example 1A. Approximate Hilbert transform pair, of orthonormal 
wavelet bases with N = 12, A" = 4,L = 2, and raid-phase spectral 
factorization. 



TABLE II 

Coefficients for Example 1 A and Example 2 



Example 1A: N = 12,tf = 4, L = 2 



*o(") 



-0.00178533012604 
0.01335887348208 
0.03609074349777 

-0.03472219035063 
0.04152506151211 
0.56035836869365 
0.77458616704024 
0.22752075128211 

-0.16040926912642 

-0.06169425120853 
0.01709940838890 
0.00228522928787 



9o{n) 



-0.00035706602521 
-0.00018475350525 
0.03259148575321 
0.01344990160212 
-0.05846672525596 
0.27464307660380 
0.77956622415105 
0.54097378940769 
-0.04031500786642 
-0.13320137936114 
-0.00591212957013 
0.01142614643933 
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Example 2: jV = 12, K = Z, L = 3 



Ao(n) 



-0.01558262447444 
-0.04943224834056 
0.21675411650608 
0.74585008428240 
0.61333711629573 
-0.01550639700556 
-0.12705042512607 
0.03236969097201 
0.01970114139115 
-0.00619091208250 
-0.00005254340590 
0.00001656336077 



9o(n) 



-0.00222608921063 
-0.04267917713309 
0.02482915969003 
0.49827824107483 
0.79972651593977 
0.28678636149680 
-0.15642754715911 
-0.03318989637197 
0.04342764217365 
-0.00220469140539 
-0.00222290024716 
0.00011594352537 



is closer to zero for negative frequencies. At the same time, the 
Sobolev exponent decreases to 1.736. This is to be expected, 
as we have reduced the number of vanishing moments arid at 
the same time increased the degree of approximation for the 
half-sample delay. The coefficients are given in Table II. 



Fig. 3. Example 2. Approximate Hilbert transform pair of orthonormal 
wavelet bases with N = 12, K = 3,L = 3„ and mid-phase spectral 
factorization. ' • ■ ' . . , 



A. Directional 2-D Wavelets 

One of the important applications of a Hilbert pair of wavelet 
bases is the implementation of directional 2-D (overcomplete) 
wavelet transforms, as illustrated in [7]. The directional 
Vavelets are obtained by first defining a 2-D separable wavelet 
basis via 

In addition, define ip 9ii similarly. Then, the six wavelets defined 
by 

M*, V) = tl>h,i{x, y) + V P> i(x, y) (9) 
^»+3(z, y) = y) - ip g>i (x, y ) (io) 
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Wn>*FHAEE CASE {EXAMPLE 1 A) 




Fig. 4. Two-dimensional wavelets generated by an approximate Hilbert 
transform pair of 1-D wavelets, using a mid-phase spectral factor (Example 
1 A) and a minimum-phase spectral factor (Example IB), respectively. 

for 1 < i < 3 are directional, as illustrated in Fig. 4. A wavelet 
transform based on these six wavelets can be implemented by 
taking the sum and difference of two separable 2-D DWTs. The 
resulting directional wavelet transform is two-times redundant. 
The inverse requires taking the sum and difference, dividing by 2, 
and the separable inverse DWTs. The four-times redundant DWT 
presented by Kingsbury is both directional and complex [8]. 

We note that the type of spectral factorization performed in 
the design procedure described above influences the quality of 
the directionality of the 2-D wavelets. For example, the wavelets 
of example 1 A, which were constructed with a mid-phase spec- 
tral factor, lead to the six 2-D wavelets shown in the top panel 
of Fig. 4. On the other hand, the wavelets of example IB, which 
were constructed with a minimum-phase spectral factor, lead to 
the six 2-D wavelets shown in the lower panel of Fig. 4. In this 
case, the directionality is not as clean as in the first case. Instead, 
some curvature is present. The figure shows that the mid-phase 
spectral factor can be preferable to the minimum-phase one for 
the implementation of directional wavelet transforms. 

Evidently, better directional selectivity is obtained when the 
two wavelets have approximate (or exact) linear phase in ad- 
dition to forming an approximate Hilbert transform pair. The 
procedure described above imposes a condition on the phase of 
of ^h(^)/^g(<^)y but it does not impose any condition of the 
phases of ^^(w) and vE^u;) individually. To ensure a high di- 
rectional selectivity, one could impose directly that H 0 (uj) and 
Go (a;) have approximately linear phase rather than relying on 
the near linear-phase of a mid-phase spectral factor. For ex- 
ample, consider the system described in [9]. In [9], the filters 
ho(n) and go(n) are related through a reversal 

9o(n) = h 0 (N - 1 - n) 



or equivalently 

Combining with (4), we get the condition 

i/oM=^" iw((A, '" 1) "' ) 

or 

H *( u ) = e -;u,((N-i)-i) 

Therefore 

H 0 (u>) = li/oHle-H^-i). (11) 

In this case, H 0 (u>) has linear phase, and the delay is offset by 
one quarter of a sample from the center of symmetry of a length 
N filter h 0 (n). As a CQF filter can not have exact linear phase, 
the filters given in [9] approximately satisfy (11). The design of 
orthonormal wavelets with approximate linear phase with non- 
integer delay has also been described in [13], [19], and [25]. An 
alternative, simple, way to obtain an approximate Hilbert trans- 
form pair of (biorthogonal) wavelets with approximate linear 
phase is to modify the spectral factorization approach, as de- 
scribed in the next section. 

III. Biorthogonal Solutions 

The quality of the directional 2-D wavelets derived from a 
pair of 1-D wavelets appears to depend on the wavelets having 
approximately linear phase, in addition to forming an approxi- 
mate Hilbert transform pair. If biorthogonal wavelets are accept- 
able, then the procedure given in Section II can be modified to 
yield wavelet bases for which both ^(w) and ^ 9 {uj) have ap- 
proximately linear phase. They can then be used to generate 2-D 
wavelet bases with improved directional selectivity. To generate 
wavelets with approximate linear phase, we can follow the ap- 
proach of [4] for the design of symmetric biorthogonal wavelets, 
in which the spectral factorization of a halfband filter is replaced 
by its factorization into two linear-phase, but different, filters. In 
the biorthogonal case, we denote the dual scaling functions and 
wavelets by ^(t), ^(0? anc * i>g(t)- As we nave a P a * r °f 
biorthogonal wavelet bases, we have eight filters including the 
dual filters, corresponding to the filterbanks illustrated in Fig. 5. 

The dual scaling function <j>h(t) and wavelet t/>h(£) are given 
by the equations 

$ h {t) = y/2Y,ho{n)$ H {2t-n) 

n 

4> h (t) = V2^2 £i(n)£ fc (2t - n). 

n 

The dual scaling function 4> g and wavelet rf> g are similarly de- 
fined. 

The goal will be to design the filters so that both the pri- 
mary (synthesis) and dual (analysis) wavelets form approximate 
Hilbert transform pairs 

= and = 



MimmrTii^wsG case <g*/uriL£ iq> 
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Fig. 5. Filterbank structures for the bibrthogonal case. 



If we define the product filters as 

Ph(n) := h 0 (n) * h 0 (n) 
■ ' Pg( n ) 9o{ n ) * 9o(n) 

then the biorthogonality conditions can be written as 

p h (2n + n 0 ) = 8(n) 
Pg(2n + n 0 )=8(n). 



(12) 
(13) 



That is to say, p/ 4 ,p<, must both be halfband. With out loss of 
generality, we can assume that n 0 is an odd integer. In that case, 
the highpass filters are given by the following expressions [4], 
[24]: 

fti(n) := (-l) n h 0 {n) 
h(n) := -(-l) n h 0 (n) 
9i(n):={~i) n g 0 (n) . '\ 

h(n) := -(-l)^o(n). 

We will look for solutions of the form 

^o(^) = f(n) * d(n) 
h 0 {n) =z f(n) * d(L - n) 
9o{n) = f(n) * d(Z, - n) 
<7o( n ) = /(n) * d(n). 

The problem is to find f (n) and /(n) such that we have have 
the following. 

1 ) The biorthogonality conditions (12) and ( 1 3) are satisfied. 

2) The wavelets have K vanishing moments (K vanishing 
moments for the dual wavelets). 

3) The wavelets form an approximate Hilbert transform pair. 

To obtain the half sample delay needed to ensure the approx- 
imate Hilbert property, we choose d(n) to come from the flat 
delay allpass filter, as before: 

A . x z- L D(l/z) U9 
A ( z ) = rv Y w z around z = l. 

To ensure the vanishing moments properties, we take F(z) and 
F(z) to be of the form 

F(z) = Q(z)(l + z-i) K 



where K + K is odd. K denotes the number of vanishing mo- 
ments of the primary (synthesis) wavelets, and K denotes the 
number of vanishing moments of the dual (analysis) wavelets. 
The product filters pn and p g are then given by 

P{z):=P h (z)=P 9 (z) (14) 
= Q«Q(*)(1 +z- 1 ) K+k D(*)D{l/z)z- L . (15) 

To obtain the required halfband property, we find a symmetric 
odd-length sequence r(n) so that 

R(z)(l + z- 1 ) K + k D(z)D(l/z)z- L 

is halfband. The symmetric sequence r(n) can be obtained by 
solving a linear system. of equations as in Section EL We can 
then obtain Q(z) and Q(z) by factoring R(z) 



R(z)=Q(z)Q(z) 



(16) 



where both q(n) and q(n) are symmetric. As q(n) and q(n) 
are symmetric, so are f(n) and /(n). It follows that h 0 (n) and 
9o(n) are related by a reversal 



and similarly 



9o(n) = h 0 (N - 1 - n) 



9o(n) =h 0 (N - 1 -n). 



(17) 



(18) 



Therefore, h 0 ,ho,g 0 and # 0 have approximately linear phase 
(because d(n) does) in addition to satisfying (4) approximately. 
Note that the symmetric factorization (16) is not unique— many 
solutions are available. In addition, note that the sequences f(n) 
and f(n) do not heed to have the same length. 
, Example 3: With K = K = 4 and L = 2, we can take the 
synthesis filters ft 0 (n) and g 0 (n) to be of length N = 13. The 
analysis filters h 0 (n) and g 0 (n) will then be of length 1 1. Fig. 6 
illustrates this solution obtained from a symmetric factorization. 
The plots of | * h (w) + j $3 (w) | and likewise \^ h (co) + j V g (w) \ 
show that they approximate zero for w < 0, as expected. The 
coefficients h 0 (n) and h 0 (n) are given in Table HI. The coef- 
ficients g 0 (n) and go(n) are given by their reversed versions, 
as in (17) and (18). The Sobolev exponent of $ h (t} and-^ p (t) 
is 1.731, whereas the Sobolev exponent of ip h (t) and tyJt) is 
2.232. 

Fig. 7 illustrates the analysis and synthesis directional 2-D 
wavelets derived from the 1-D wavelets using (9) and (10). 

IV. IIR Solutions 

The spectral factorization approach can also be used to con- 
struct orthogonal wavelet bases based on recursive infinite im- 
puilse response (HR) digital filters, where H 0 (z) is a rational 
function of z. Wavelets based on rational scaling filters have 
been discussed, for example, in [6], [15], [17], and [21]. UR 
filters often require lower computational complexity than finite 
impulse response (FIR) filters. Analogous to the approach given 
in [6], but with the flat delay filter included, approximate Hilbert 
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TABLE III 

COEFFICIENTS FOR EXAMPLE 3 

Example 3: AT = 13, K = X = 4, L = 2 



Mn) 



M") 



0.01339704408541 
-0.00678912656592 
-0.15762026783144 
-0.04891448800028 
0.64995392445932 
0.93376892975667 
0.27109270647926 
-0.19103597922739 
-0.07239603482308 
0.02007744522348 
0.00267940881708 
0 



-0.00030453648331 

-0.00015432782903 
0.02776733337848 
0.01114383888330 

-0.04715783038404 
0.23730846518842 
0.65840893759976 
0.47625050860818 
0.03941149716348 

-0.02885152385159 
0.03050406232874 
0.01140982018728 

-0.00152268241655 
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Fig. 6. Example 3: Approximate Hilbert transform pair of biorthogonal 
wavelet bases with N = 13, A* = K = 4, L = 2. 



transform pairs of ILR wavelets can be obtained with the fol- 
lowing form: 



Fig. 7. Two-dimensional wavelets generated by an approximate Hilbert 
transform pair of biorthogonal 1-D wavelets (Example 3). 



The product filter is given by 

P(z) := H 0 {z)H 0 (l/z) = Go(*)G 0 (l/z) 
(z + 2 + z- 1 ) A 'D(z)Z>(l/z) 



G 0 (z) 



(1 + 


z" 1 ) A 'Z?(z) 




C(*») 


(1- 


z- l ) K D(-l/g)z- 




C(»») 


(1 + 






C(z») 


(1- 







Defining 



C(z 2 )C(l/z 2 ) 



VW := (z + 2 + z- 1 ) A D(z)D(l/z) 



the orthogonality condition P(z) + F(— z) = 2 can be written 
as 



(19) 
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Fig. 8. Example 4. Approximate Hilbert transform pair of orthonormal IIR 
wavelet bases with K = 2, iL = 2. " . 



C(^) can be found by spectral factorization; note that the 
left-hand side of (19) is a function of z 2 . A stable filter is 
obtained using the minimum-phase spectral factor in ( 1 9). 

Example 4: With two vanishing moments (K = 2) and L = 
2, we obtain the following stable causal transfer functions 



#o(*) = 
G 0 (z) = 



1 + fe" 1 + 5:2*-? + 2.4z- 3 '+0.2^- 4 
. 6.6495 H- 2.3714z- 2 -f0.0301z- 4 
0.2 + 2.4*" 1 + 5.2z~ 2 ,+ 4?- 3 + z~ 4 
6.6495 ■+ 2.3714^- 2 + 0.0301z~ 4 * 



Fig. 8 illustrates the two' wavelets iph{t),i> g (t) and the magni- 
tude of *fc(w)+j*^(ci;). 

V. Comparison to Kingsbury's Filters 

It is interesting to compare the filters obtained here with 
those presented by Kingsbury in [7]-[10], where the com- 
plex dual-tree DWT is introduced. In [7] and [8], Kingsbury 
presents linear-phase biorthogonal solutions where one scaling 
function is symmetric about an integer, whereas the other 
scaling function is symmetric about an integer plus 0.5. In 
-this case, the (approximate) half-sample delay [see (4)] can 
be achieved by making the frequency responses magnitudes 
(approximately) equal -\H 0 (u>)\ » |G 0 (o;)|. The family of 
filter banks introduced in [9] and [10] are orthogonal but with 
the property that one scaling function is the time-reversed 
version of the other. Then, the frequency responses magnitudes 
are equal, and it remains to design the phase response of the 
filter to achieve (approximately) the half-sample shift [see 
(11)]. On the other hand, the approach described in this paper 
does not impose any condition on the phase-response of the 
filters. The orthogonal solutions given in this paper do not have 
approximately linear phase (for short filters the mid-phase 



spectral-factorization does not always provide a solution with 
a nearly linear phase response). However, the biorthogonal 
solutions given in this paper have approximately linear phase. 
As a result, the orthogonal solutions by Kingsbury are more 
nearly linear phase than the orthogonal solutions given here; 
likewise, the linear-phase biorthogonal solutions by Kingsbury 
are more nearly orthogonal than the approximately linear-phase 
biorthogonal solutions given here. 

In addition, the design procedure developed in this paper 
and the ones described by Kingsbury are quite different. 
Kingsbury's design procedures are based on minimizing the 
aliasing of a filterbank after it is iterated a fixed number stages. 
The minimization is performed over a parameterized space of 
perfect reconstruction filterbanks using iterative optimization 
algorithms. On the other hand, the approach taken in this paper 
is based on the limit functions [defined by (2) and (3)] and on 
vanishing moments. As a result, the limit functions Vv* (*) and 
i> g (t) are closer to forming a Hilbert pair for the filters given 
in this paper than those presented by Kingsbury, whereas the 
solutions by Kingsbury do this better for the first several stages. 

VI. Conclusion 

This paper presents simple procedures for the design of pairs 
of wavelet bases where the two wavelets form an approximate 
Hilbert transform pair. The approach proposed here is analo- 
gous to the Daubechies construction of compactly supported 
wavelets with vanishing moments but^ where the approximate 
Hilbert transform relation is added by way of incorporating a 
flat delay filter. 

The approach is based on a characterization of Hilbert trans- 
form pairs of wavelets bases given in [18]. The formulation of 
the problem, using a flat delay allpass filter, makes it possible to • 
employ the spectral factorization design method as introduced 
in [22] for the design of CQF filters. Note that even though an 
allpass filter arises in the problem formulation, the filters we 
obtain are FIR (HR solutions are also available, as described in 
Section IV). 

Given an allpass filter, the proposed design method produces 
short filters with a specified number of vanishing wavelet mo- 
ments. Although a flat delay filter was used here, any other all- 
pass filters that approximate a delay of a half sample could be 
used instead. The degree of the allpass filter controls the degree 
to which the half-sample offset property is satisfied. A Matlab 
program for the procedure is given, and examples are also given 
to illustrate the results. Programs are available on the Web at 
http://taco.poly.edu/selesi. 
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